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Abstract 

A  molecular  dynamics  simulation  of  liquid  nitromethane  has  been  carried  out,  in  which  a 
density  functional  (DF)  procedure  was  used  to  compute  the  properties  needed  to  establish  the 
intermolecular  force  field.  The  latter  was  subsequently  updated  by  further  DF  calculations.  The 
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first  and  second  shells  of  neighbors  are  found  to  be  at  radii  of  about  6  A  and  1 1  A  from  a  given 
nitromethane  molecule.  OH  partial  pair  correlation  functions  are  determined,  and  indicate  the 
presence  of  C-H . O  hydrogen  bonding. 
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Nitxomethane,  a  prototypical  energetic  compound,  is  a  liquid  at  room  temperature,  with  a 
boiling  point  of  101.2°C  at  one  atmosphere  pressure.^  The  polar  nature  of  the  molecule  is 
reflected  in  its  dipole  moment,  3.46  D,  and  in  the  dielectric  constant  of  the  liquid,  35.87  at  30°C.^ 
Quite  significant  intermolecular  interactions  can  therefore  be  anticipated  in  the  liquid  phase,  perhaps 

including  C-H . O  hydrogen  bonding.  (It  is  now  well-established  that  carbon  can  be  a  hydrogen 

bond  donor.  )  The  elucidation  of  such  interactions  would  be  important  both  in  furthering  the 
understanding  of  the  structures  of  polar  liquids  in  general  and  also  in  helping  to  interpret 
experimental  studies  of  nitromethane  in  particular.^'^^ 

We  have  accordingly  undertaken  a  combined  density  functional  (DF)/molecular  dynamics 
(MD)  study  of  liquid  nitromethane  at  300°  K  and  a  density  of  1 . 139  g/cm^.  We  use  a  DF 
procedure,  applied  to  a  single  molecule,  to  compute  the  properties  needed  to  establish  the 
intermolecular  force  field  and  then  begin  an  MD  simulation.  This  is  stopped  occasionally  and  a 
new  DF-updated  force  field  is  obtained,  using  now  a  typical  configuration  of  several  molecules. 
This  continues  until  self-consistency  is  attained  in  the  force  field,  which  is  then  appropriate  to  the 
specific  conditions  of  temperature  and  density.  Since  no  experimental  data  are  required  as  input, 
this  approach  can  be  applied  as  well  to  compounds  which  have  not  yet  been  synthesized  or  for 
which  measurements  have  not  been  carried  out. 

The  molecular  properties  needed  to  create  the  initial  force  field  are  the  optimized  geometry, 
force  constants  and  atomic  charges.  Once  the  MD  simulation  is  underway,  the  periodic  DF 
calculations  on  the  nitromethane  clusters  are  single  point  runs  to  determine  new  charges;  the 
geometries  used  are  those  that  have  evolved  during  the  previous  MD  simulation.  Our  approach  can 
be  viewed  as  an  improvement  over  MD  techniques  that  use  a  non- varying  force  field  (which  may 
be  derived  from  ab  initio,  semi-empirical  or  experimental  data)  and  in  the  direction  of  the  procedure 
of  Car  and  Parrinello,  which  directly  combines  the  DF  and  MD  equations  so  that  no  force  field 
need  be  explicitly  created.  While  the  latter  approach  is  certainly  a  very  desirable  one,  it  is  limited 
at  present  to  relatively  small  systems. 

The  density  functional  calculations  were  carried  out  with  Gaussian  92/DFT,  using  the 
Becke  exchange^^  and  Lee-Yang-Parr  correlation^^  functionals  and  the  6-3 IG**  basis  set.  Atomic 
charges  were  obtained  by  fitting  them  to  the  computed  molecular  electrostatic  potential  while 
simultaneously  preserving  the  dipole  moment.'^  The  resulting  charges  were:  N:0.660;  O:  -0.377; 
C:  -0.266;  H:  0.120.  The  CHARMM  program  was  used  for  the  molecular  dynamics 
simulations^^  and  analysis  programs  were  developed  in  situ.  216  nitromethane  molecules  were 
distributed  uniformly  and  with  random  orientations  in  a  cubic  box  with  dimension  26.79  A,  chosen 
to  reproduce  the  experimental  density.  ^  Periodic  boundary  conditions  were  used.  The  dynamics 
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were  performed  with  the  Verlet  algorithm. 


Before  beginning  the  simulation,  a  classical  optimization  of  the  orientations  and  geometries 

of  the  molecules  was  carried  out.  The  resulting  average  geometry  was  very  similar  to  that 

produced  by  the  DF  calculation,  the  largest  difference  being  a  shortening  of  the  C-N  bond  by  0.02 
0 

A.  In  the  early  stages  of  the  simulation,  the  temperature  of  the  system  was  increased  by  10°K 
every  1  ps  (2000  steps);  later  this  became  10°  K  every  2.5  ps  (5000  steps).  Once  the  final 
temperature  of  300°K  was  reached,  the  system  was  allowed  to  equilibrate  for  20  ps.  During  the 
final  portion  of  the  equilibration,  the  potential  energy  and  the  temperature  showed  only  minor 
fluctuations. 

Some  time  after  equilibration  had  been  finished,  at  52  ps  after  the  beginning  of  the 
simulation,  a  cluster  of  nine  nitromethane  molecules  was  selected  for  the  DF  update  of  the  force 
field.  This  cluster  consisted  of  a  central  molecule  and  its  first  shell  of  neighbors.  The  new 
charges  found  for  the  central  molecule  (and  adopted  for  all  216)  were;  N:  0.717;  O:  -0.406;  C: 
-0.33 1 ;  H:  0. 142.  A  subsequent  DF  analysis  after  20  ps  (40,000  steps)  showed  no  significant 
further  changes  in  atomic  charges.  Once  the  system  had  adjusted  to  the  modified  charges,  desired 
properties  were  determined  by  taking  samples  at  various  regular  intervals,  using  a  step  time  of  0.5 
fs. 

Figure  1  shows  the  CC,  CN  and  NN  pair  correlation  functions  as  determined  from  the 
production  runs.  These  show  that,  for  a  given  nitromethane  molecule,  the  first  shell  of  neighbors 
is  within  a  radius  of  about  6  A,  and  the  second  shell  is  within  about  1 1  A.  Of  particular  interest  are 
the  OH  correlations,  shown  in  Figure  2.  For  any  pair  of  nitromethane  molecules,  there  are  12 
possible  0-H  interactions.  The  total  OH  pair  correlation  function  is  the  average  of  all  of  these, 
resulting  in  the  essentially  unstructured  curve  near  the  bottom  of  Figure  2.  We  were  able  to  resolve 
this  into  its  12  components  (not  readily  achievable  experimentally)  that  correspond  to  the  12 
possible  0-H  intermolecular  distances  between  two  nitromethane  molecules.  The  first  peak  of  the 
first  OH  partial  pair  correlation  function  occurs  at  about  2.3  A;  this  is  clearly  indicative  of  C-H-0 
hydrogen  bonding  in  liquid  nitromethane,^.^  a  point  which  has  been  the  subject  of  speculation  in 
the  past. 

The  self-consistent  force  field  procedure  that  we  have  applied  in  this  study  provides  a  more 
realistic  and  yet  computationally  efficient  means  of  treating  large  systems.  In  addition,  the  use  of 
partial  pair  correlation  functions  permits  a  much  more  detailed  analysis  of  the  structure  of  fluids 
than  is  possible  with  the  total  pair  correlation  function. 
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Captions  for  Figures 

Figure  1.  Calculated  CC,  NN  and  CN  pmr  correlation  functions  for  liquid  nitromethane  at  300  °K 
and  a  density  of  1 . 139  g/cm^ 

Figure  2.  Calculated  OH  total  (thick  curve)  and  partial  pair  correlation  functions  for  liquid 

nitromethane  at  300  °K  and  a  density  of  1.139  g/cm^.  The  12  partial  pair  correlation 
functions  correspond  to  the  possible  O — H  intermolecular  distances  between  two 
nitromethane  molecules. 
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